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ABSTRACT 


The purpose of this study is to determine the magnitude of 
error introduced due to wing vibration when measuring atmospheric 
turbulence with a wind probe mounted at the wing tip and to determine 
whether accelerometers mounted on the wing tip are needed to correct 
for this error. A spectrum analysis approach is used to determine 
the error. Estimates of the B-57 wing characteristics are used to 
simulate the airplane wing, and von Karman's cross spectrum function 
is used to simulate atmospheric turbulence. The major finding of 
the study is that wing vibration introduces large error in measured 
spectra of turbulence in the frequency's range close to the natural 
frequencies of the wing. 
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1.0 INTRODUCTION 


Progress in applying spectrum analysis in aeronautical engi- 
neering is fostering the measuring of the atmospheric turbulence 
cross spectra, especially across wing spans. Measuring the atmo- 
spheric turbulence spectrum across the wing span requires wind sensors 
mounted on the wing tips, where the sensors are exposed to error due 
to wing vibrations. Determining this error requires comparison 
between the spectrum of the wing tip velocity and spectrum of 
atmospheric turbulence. This study demonstrates a procedure to 
estimate the wing tip velocity spectrum and the error introduced due 
to wing vibration when measuring atmospheric turbulence with a wind 
probe mounted at the wing tip of a B-57 type airplane. The wing 
characteristics of a B-57 are used for modeling the airplane wing 
because NASA is currently planning to use this airplane to acquire 
atmospheric gust gradient data across the wing span. The purpose of 
this study is to determine whether the error introduced in measuring 
atmospheric turbulence is large enough to justify mounting accelero- 
meters on the wing tip to correct for the error introduced into the 
measured turbulence data due to the wing vibration. 

An introduction to spectrum analysis is believed useful in 
order to clarify the general solution technique and introduce the 
major terms of spectrum analysis. Spectrum analysis is based upon 
transferring the system from the time domain to the frequency domain. 

The form of the response of a linear system in the time 
domain to a single arbitrary input p(t) is given by Duhamel's integral: 
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w(t) = 


CO 


P(t) h(t - T)dr 


(1) 


J *00 


where h(t) is the response of the system to a unit impulse. If p(t) 
is a random input, evaluating the integral becomes meaningless because 
the integral represents only one sample of a random distribution of 
responses. To overcome this difficulty, Duhamel's integral can be 
transformed from the time domain to the frequency domain for stationary 
random inputs. The Fourier transform of Duhamel's integral is: 

W(w) = P(u) Z((o) (2) 

where the frequency response function Z(w) is the Fourier transform of 
h(t) and represents the response of the system to a sinusoidal input. 

In general, the functions in Equation (2) are complex valued functions 
and contain both magnitude and phase information. In order to study 
only the magnitude of the response in the frequency domain, the power 
spectrum of the response is defined by: 

c|)^(aj) = W(co) W{o)) = P(u)) Z(oj) P(w) Z(to) = cj)p(a)) lZ(co)l^ (3) 

where <f>p is the power spectrum of the input. This equation is the 
fundamental result of spectral analysis and equates the response spec- 
trum to the product of the input spectrum and the square of the 
magnitude of the frequency responsco References [1] and [2] discuss the 
derivation of the theory and contain many examples of the application 
of spectrum analysis. 
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The input-output relationship of Equation (3) can be used to 
outline the general solution technique of spectrum analysis. First, 
the frequency response function, Z, is determined via the "equation 
of motion"; second, the input spectrum, <|)p, is defined; and, third, 
the product in the right-hand side of Equation (3) is calculated to 
obtain the output spectrum. 

The complicated system involved with wing vibrations requires 
simplification in order to become amenable to solution. The major 
simplification in determining the frequency response function is 
putting restraints on the aircraft. This study will restrain the 
aircraft to only vertical motion of the center of gravity and vertical 
wing bending. The coordinate system used in this study is shown 
in Figure 1, where the spanwise coordinate, y, is the independent 
variable and the vertical deflection, w, is the dependent variable. 

The aerodynamic forces appearing in the equation of motion are those 
calculated from strip theory and assume the wing is a flat plate 
having stiffness similar to a B-57 wing. The gust is considered to 
include only vertical wind variations. The atmospheric spectrum is 
assumed to be stationary, homogeneous and isotropic; but turbulence 
varies across the wing span. For calculation purposes, the von Karman 
cross spectral function as defined by Houbolt and Sen [3] is assumed. 
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2.0 ANALYSIS 


2.1 Spectrum Analysis 

The spectrum equation for a single stationary random input, 
derived in the introduction, is not sufficient for calculating the 
spectrum of the wing tip velocity for an aircraft with large wing span. 
Assuming that the aircraft experiences a single input or gust would be 
tantamount to assuming that the gust field is uniform across the span 
but random in the flight direction, as depicted by the left-hand sketch 
in Figure 2, This assumption is unrealistic for an aircraft 
with large wing span. A better assumption is that the gust is two- 
dimensional, so that it is also random across the span of the aircraft, 
as depicted in the right-hand sketch in Figure 2. The assumption of 
one-dimensional turbulence may lead to an underestimation of the 
response of the aircraft. Houbolt [4] has shown the the root bending 
moment spectrum determined by two-dimensional turbulence is signifi- 
cantly higher than the root bending moment spectrum for one-dimensional 
turbulence. The ratio 2£/L, where 2z is the span of the aircraft and 
L is the turbulence length scale, detennines the validity of the assump- 
tion of uniform spanwise turbulence; for 2z/L approaching unity the 
assumption of spanwise uniformity becomes invalid. The B-57 has a wing 
span of 2sl = 66 ft, and Houbolt [5] recommends a length scale of atmo- 
spheric turbulence of L = 300 ft, in which case the ratio 2£/L = 0.22 
is large enough to warrant a two-dimensional spectrum analysis. Near 
the ground L becomes smaller making the above argument even stronger. 
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FIGURE 2 ONE-DIMENSIONAL AND TWO-DIMENSIONAL TURBULENCE ENCOUNTER [4] 

Assuming a linear system with a continuum of random stationary 
inputs, the response in the form of Duhamel's integral is: 

, £ - oo 

w(t) = p(y.T) h(y,t - t) dr dy (4) 

■ -£ •' -°o 

The Fourier transform of the above equation is: 

a 

W{o)) = P(y.o)) Z(y,(o) dy (5) 

■ -Z 

where the frequency response function Z(y,to) is a function of both 
frequency w and input location y. This study will define Z(y,co) as 
the velocity of the right wing tip due to a sinusoidal gust of 
frequency u located at y along the wing and w(t) as the velocity of 
the right wing tip due to gust excitation along the entire wing. 

The output power spectrum as defined by Equation (3) is 

j- £ /■ 5, 

(j)„((o) = W(w) ¥(w) = <|)p(yi .y2.u) Z(y^,(j) ZCyg.u) dy^ dy2 

• -Z J-£ 
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( 6 ) 
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The cross spectrum is defined by: 

4>p(yi>y2»^) ~ P(yi»m) P (y2 »u) (7) 

Then it is true for all cross spectra that 

4>p(yi»y2»^) ~ PCyisfij) P(y2 

~ P (y"! >w) P(y 2 >w) 

= ?p(y 2 .y-|>w) ( 8 ) 

Because the domain of integration is symmetric, an appropriate choice 
of the limits of integration can be made so that the integrand can 
be written as: 

'f)p(yi .y2>w) Z(ypu) Z(y2,w) + (f>p(y2. y-|.w) Z(y2.to) Z (y-, ,to) 

.= <(>p(yi >y2>“) Z(y^ ,w) Z(y2>(jo) + ({jpCy-] ,y2>w) Z(y2,to) Z(y^ ,w) 

= <J>p{y-| .y2>w) Z(y^,0)) Z(y2,w) + <j>p(yi ,y2,o)) Z(y^ ,0)) Z{y^,oi) 

= 2Re[(|)p(y^ ,y2,cj) Z(y^ ,o)) Z(y2,w)] (9) 

Therefore the symmetry involved in the domain of the integration and 
the products of conjugate pairs permits the integral to be rewritten as: 

• i 

<})p(yi .y2»“) Z(y^,0)) 7(y2,u)) dy^dy2 

y* 

( 10 ) 

7 


= 2 Re 


lim 

y*-^yc 


L 



This equation is completely general and requires only that the inputs 
be stationary. It is reasonable to expect that the spectrum of the 
input (in particular, turbulent input) is isotropic; in other words, 
the spectrum is a function of the separation distance [meaning 
^p^^l’^2^ “ '^p^l^l " ^ 2 !^^ independent of direction. This 

assumption can further simplify the integral by making the variable 

substitutions s = y^ ~ ^2 ^ ^ ^2 interchanging the order of 

integration so that: 


lim 

r 25, 


• 5,-S 

4)y^(w) = s*->0 

4) (s,o)) 2Re 


7(y,w) Z(y + s,(jo) dy 

• 

s* 

. 

-5, 


This equation is the final form of the spectrum equation that 
will be used in this study. Note that (j)p(s,o)) is the cross spectrum of 
the inputs, which within this study is the cross spectrum of atmospheric 
turbulence. The cross spectrum of turbulence represents the contribu- 
tion of the frequency w to the covariance of the vertical velocity at 
two points along the wing span separated by a distance s. 

2.2 Mechanical Analysis 

Consider a nonuniform wing (chord, beam stiffness, and mass 
varying along the span) free to move and bend vertically and subjected 
to a sinusoidal vertical gust at a spanwise element of width A centered 
at y = y*. The balance of forces then requires that 

Fs = -Fi - F^ + Fg 6(y,y*) (12) 
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where is the force due to beam stiffness, Fj is the inertial force, 
and F|^ and Fg are aerodynamic forces due to translatory wing motion and 
due to vertical gust, respectively. The function 6(y,y*) selects the 
portion of the wing which is subjected to the gust and is zero every- 
where except between y* ± A/2, where it has the value of unity. The 
theory of strength of material states: 


F 


S 



El 



(13) 


where El is the beam bending stiffness and w is the vertical deflection 
of the wing. Newton's law gives: 


F 


I 


= mw 


(14) 


Theodorsen in his famous NACA report [6] showed that the force due to 
wing motion for a two-dimensional wing is: 


F^ = Trpb^w + 2frpVb C(k) w 


(15) 


where V is the mean flight speed, p is the density of the flight medium, 
and b is the semi -chord of the wing. The Theodorsen function, C(k), is 
a function of the reduced frequency (k = wb/V) of the motion. Both 
References [2] and [7] develop and apply the Theodorsen function. 

Defined in terms of the Bessel function, the Theodorsen function is: 

. ■1i(Ji^Yo)^Y^(V, -Jo)-i(Y,Y, 

(J, t Y„)2 + (Y, - 
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Note that the argument of the functions is k, the reduced frequency of 
the motion. The lift due to a gust acting on a two-dimensional wing 
with semi -chord b as reported in [2] is: 


Fg = 2iTpV^b 


K(k) 


(17) 


where K(k) is the Kussner function. It too is dependent on reduced 
frequency of the gust. In terms of the Theodorsen function and Bessel 
functions, the Kussner function is defined: 


K = C(Jq - iJ^) + iJ^ 


(18) 


The forces may be added to yield the partial differential equation: 


3^w 


El 

9y" I 9y"J 


= -mw - Trpb^w - 2iTpVb C(k) w + 2irpV^b K(k) ^ 5(y.y*) 


(19) 


The boundary conditions are; 


w"(£,t) = w"'(£,t) = w"(-£,t) = w’"(-£,t) = 0 (20) 

meaning no shear or moment at the wing tips, while the time dependence 
will be taken to be periodic, as discussed later. 

Assume the vertical deflection w(y,t) may be expanded in its 
natural modes, as: 

00 

w(y,t) = X) n,-(t) (j).(y) (21) 

i=l ^ ^ 
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Substitution of the solution expanded in its natural modes into the 
governing Equation (19) yields, after some algebra, an equation for the 
spatial dependence (the free vibration equation) and another for the 
time dependence: 


9y 


El 


ay 


2 'f’i 




for 


-z < y < a 


( 22 ) 


where is the natural frequency of the ith mode. The boundary condi- 
tions of the differential equations are: 


<},."■(£) = "■{-«-) = = 0 (23) 


which corresponds to the shear and bending moment being zero at the ends, 
as is the case for free ends. 

The natural modes are orthogonal [2] because of the choice of 
boundary conditions so that: 


Z 

(Ei<p.")'’ cfjdy = 

• -Z 


Z 

m<f)^-<|)jdy = 

• -Z 


2 

for i = j 

- 

0 for i j 

(24) 


Then by using the orthogonality property of the modes, the 
differential Equation (19) becomes a system of linear differential equa- 
tions in terms of the natural modes: 


Vn^^n ' -Vn ‘ ''n/j ‘ ^ 

J I J * 


+ 2TTpV^bp K(k) ({,.(y*) a(y*)A 


(25) 
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Where 




a 


and 


B . = 
nj 


z 

a^n<Djdy 

i-z 


are coefficients of the aerodynamic cross terms, the weighting function 
of the integrals a(y) = b(y)/bp is the semi-chord distribution. We 
now introduce the variables: 


k = 


OjbR 

V 



(27) 


and assume sinusoidal output to sinusoidal input; in other words, we 

T 1^5 — 1 ks 

assume the variables are of the form u = ue and = n^e . After 

2 1 ks 

substitution of the above variables and division by upV Se , the 
system of linear differential equations becomes a system of linear 
algebraic equations: 


= k^y„5n + E + 2ik C(k) £ 


N 


N 


where n 




n 


An-! 

nj 


b„A . 
R nj 


^nj 


b„B . 
R nj 


2b. K(k) 


t|)„(y*) a(y*) 


(28) 
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M 

_ n 

Trpbj^S 



This is the final form of the equations which will be used in 
the program. Many of the constants of the equations are calculated by 
integrating products of the natural modes by different weighting func- 
tions, and they can be determined as soon as the free vibration problem 
has been solved for the different modes. The solution = ^^/bj^(V/LT) 
represents the amplitude of the modal response of the deflection of the 
wing to sinusoidal gust with unit change in angle of attack normalized 
by the reference semi -chord. The solution, may be 

differentiated with respect to time to give the amplitude of the wing 
tip velocity of the modal response: 


d_ 

dt 


""n V _ , % 


( 29 ) 


For the particular problem of this study, the frequency response 
function is defined as the velocity of the right wing tip due to a 
sinusoidal gust located at y* along the wing. The frequency response 
function is then: 


Z(y*,u) 


N 


= io) E 
j=l 


nj(y*,w) 

F 


( 30 ) 


This is the definition used in the spectrum analysis part of the program. 
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3.0 NUMERICAL PROCEDURE 


The numerical procedure is basically divided into three sub- 
routines. The first subroutine solves the free vibration problem and 
determines the natural bending modes and frequencies of the wing. The 
second subroutine solves the forced vibration problem and determines 
the frequency response function of the wing to a sinusoidal gust at a 
point along the span. Finally, the third subroutine performs the arith- 
metic in the power spectrum equation and determines the wing tip velocity 
spectrum. 

3.1 Free Vibration: Structural Character of the Ming 

The eigenvalue problem solved by this program is for a wing 
vibrating in only bending modes in free space (meaning free ends, i.e., 
moment and shear are zero at the wing tips). 

The differential equation is: 

(EI(J)")" = mu^(f> for -Z < y < Z (31) 

with boundary conditions: 

rOi) = <f)"(il) = ((,•■•(-£) = f’(-£) = 0 (32) 

where 4> is the natural mode or eigenfunction and u is the natural fre- 
quency or eigenvalue. 

The coefficients of this equation are El, the bending stiffness 
of the wings, and m, the mass per unit length, and both are a function 
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of y. Both El and m are defined syrnmetrically for most wings, and this 
implies that the solution of the differential equation should be either 
symmetric or antisymmetric. The program takes advantage of this so that 
the range of integration is cut in half (i.e., the problem needs to be 
solved only for 0 < y < 2.). The boundary conditions must then be 
restated at the origin (or the midspan of the airplane) to be those for 
a symmetric or antisymmetric function dependent upon whether an even or 
odd mode is being investigated. The conditions at the origin for a 
symmetric function are: 


(|)'(0) = = 0 (33) 

while the conditions for an antisymmetric function are: 

(j)(0) = ((>"(0) = 0 (34) 

These will be the boundary conditions, along with those at the wing tip 

(y = 2-). 

If the eigenvalue of the differential equation is specified, 
the problem is reduced to a linear fourth order, two-point boundary 
value problem, and a shooting method [10] can be used to determine the 
numerical solution. The first step of the shooting method is to reduce 
the single fourth order equation to a system of four first order equa- 
tions. For Equation (31) the system becomes: 


15 



■ 4*1 ■ 


0 

1 

0 

0 “ 


■ 4*1 ■ 

4*2 


0 

0 

1 

0 


4*2 

CO 

-e- 


0 

0 

0 

1 


4>3 



mo)^ 

0 

-2EI' 

-El" 


.‘^4. 


1— » 
LiJ 

El 

El 



and boundary conditions for a symmetric mode are: 


(35) 


(36) 


If all the initial conditions were given, a Runge Kutta scheme could be 
used to determine the numerical solution, but because the solution must 
match the boundary conditions at y = £, a shooting method must be used. 

The shooting method estimates the complete set of initial conditions 
and then uses a Runge Kutta scheme to determine the solution of the 
initial value problem. The value of the solution of the initial value 
problem at y = Z is compared with the boundary conditions at y = £ for the 
two-point boundary value problem. The estimate of the complete initial 
condition is improved and the procedure is repeated; in this way the 
complete set of initial conditions is determined so that the solution 
to the initial value problem has the correct value for the boundary 
conditions at y = £. Fortunately, for a linear system the complete set 
of initial conditions can be improved to the correct initial conditions 
after one trial. This will be shown true later. Because the differen- 
tial equation is linear, it can be shown [8] that there exists a vector 

— 4 / V 

base, of the solution of Equation (35). Then each solution 
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( 37 ) 


<{>1 

‘J>p 

<p = 

4)3 

of Equation (35) for the boundary conditions indicated is a linear combi- 
nation of the base solution: 


4 _ 

4> = Z ip.C. (38) 

i=l ^ ^ 

By choosing four linearly independent initial conditions, ^^^=1 is 
guaranteed to be a complete base. A convenient choice of linearly 
independent initial conditions is: 


■ 1 


' 0 ■ 


■ 0 ■ 


■ 0 ■ 

0 


1 

; 4^3(0) = 

0 

and 4^4(0) = 

0 

0 


0 


1 


0 

_ 0 _ 


. 0 . 


. 0 . 


1 


(39) 

A base of defined by the above initial conditions and the 

differential equation is commonly called a fundamental base or the fun- 
damental solutions. The boundary conditions are used to determine the 
4 

constants {C ^. }^._1 of Equation (38). In the case of the even mode, the 
4 

{C ^. }^.^1 must satisfy the system of equations: 
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^42^^^ 

11^14(0) ’^^44(0) 

'I^13(A) ^ 23 ^^^ ^ 33 ^^^ ^ 43 ^^^ 

'^24^^^ '^34^'^^ 'f’44^®) 

If the problem was a simple two-point boundary value problem, the 

boundary conditions would be nonhomogeneous. Then the system of linear 

4 

Equations (40) would be nonhomogeneous and the {C^.}^._1 could be deter- 
mined after the fundamental solutions were determined from one 

trial of the shooting method. An eigenvalue problem by definition 

requires that the boundary conditions be homogeneous, and this leads to 

4 

a homogeneous system of linear Equations (40) for the . Linear 

algebra theory requires that the determinant of the coefficient matrix 
vanish in order for the homogeneous system to have a nontrivial solu- 
tion. This determinant: 




ro 

ro 

0 

11^32(0) 

0 

CVJ 

'1^14(0) 

^ P 24 (^) 

*34(0) 

4^44(0) 



ij; 33 (£) 

'^43 ^ 

4 ^T 4 (^) 


'l^ 34 ('^) 

^44 


4 

is commonly called the characteristic determinant, and for 

to be determined it must vanish. Note that the determinant is a function 

of (0, the unspecified parameter of the differential equation. 
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Reference [9] shows that the characteristic determinant is an analytic 
function of w and the eigenvalues of the differential equation are the 
zeros of the function D(w). When u is an eigenvalue of the differential 
equation, the characteristic determinant vanishes and the boundary con- 
ditions can be satisfied as a linear combination of the four fundamental 
solutions. Note that the solution of the system of linear equations 
for the is not unique; therefore an extra condition must be sup- 

plied. In this study we demand that the natural modes have the value of 
unity at the wing tips, which becomes our extra condition imposed on 
the . The values of the derivatives of the fundamental solutions 

appearing as elements of the characteristic determinant are obtained by 
solving the differential equation by a Runge Kutta method. Note that 
some of the terms in the determinant are already known from the defini- 
tion of the fundamental solutions. Substituting in the determinant for 
these values and simplifying yields; 


D(o)) = 


^ 43 ('^) 

^ 24 ^^^ ^ 44 ^^^ 




(42) 


The characteristic determinant is now a function of only the 
first and fourth fundamental solutions because these are the only two 
fundamental solutions that satisfy the conditions at the origin for a 
symmetric function. The program takes advantage of this and determines 
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only two fundamental solutions. Reference [9] further discusses the 
theory of linear differential equations. 

The general solution procedure can now be outlined. First, 
the eigenvalue is estimated; then the two fundamental solutions are 
determined by a Runge Kutta Fehlberg order seven scheme. The value of 
the characteristic determinant is calculated from the fundamental 
solutions. A search routine checks if the eigenvalue is bracketed 
between the current estimate and the previous estimate. In this case, 
the program is directed to a bisection routine to improve the brackets 
or continues, taking another step along the frequency line and using 
this as its next estimate of the eigenvalue. After the eigenvalue has 
been determined, the natural mode is normalized by a unit displacement 
at the right wing tip. The new mode is integrated with previous modes 
determined to calculate the aerodynamic cross terms. The program then 
steps along the frequency line for its first and second estimates of the 
next eigenvalue. Another example of this technique for solving eigen- 
value problems is found in Reference [10]. 

The calculation was tested against the uniform beam and was 
found to be very accurate. Three runs were made for different El 
distributions, and the program converged very quickly for the lower 
modes. Figures 3 and 4 show estimates of m and bj^ for the B-57 used in 
all the cases run. 

The most difficult parameter to estimate is the beam stiffness, 
El. A static analysis was used to estimate the beam stiffness, assuming 
a loading on the B-57 wing which would be used during extreme operations 
and a load factor of 10 g. Appendix B shows the details of the analysis 
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FIGURE 3 DISTRIBUTION OF MASS ACROSS THE SEMISPAN OF THE WING 



0 11 22 33 

Distance from midspan (ft) 


FIGURE 4 DISTRIBUTION OF SEMI-CHORD ACROSS THE SEMISPAN OF THE WING 
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the results of the analysis show that the beam stiffness is approximately 


El = 9 X 10^ for y < 111 I 

El = 9 X 10^ for y > |11 [ (43) 

This wing is referred to as the "standard wing" throughout the rest of 
this study. Its natural modes and frequencies are shown in Figure 5„ 
From an in-flight experiment with the B-57, NASA/Langley Research Center 
determined that the natural frequency of the first bending mode of 
the B-57 is approximately 7 Hz. In this study, a trial and error 
method was used with the program to determine the beam stiffness 
necessary to have a first bending mode at 7 Hz. The beam stiffness 
was found to be: 


El = 3 x 

10® 

for 

y 5 

in 

El = 3 X 

10® 

for 

1 V 

in 


(44) 


The wing with this beam stiffness distribution is referred to as the 
"stiff wing." For comparison purposes, a third beam stiffness with 
distribution 


El = 9 X 

10® 

for 

*< 
1 A 

m 

El = 9 X 

10® 

for 

A 1 

>> 

in 


(45) 


was run. Figure 6 shows the natural modes and frequencies of this 
"flexible wing. " 
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Normalized wing displacement 


Reduced natural frequency 


Natural frequency 
Mode (Hz) 

1 4.37 0.072 

2 12.06 0.199 

3 23.35 0.385 

4 40.34 0.666 



FIGURE 5 NATURAL MODES AND FREQUENCIES FOR THE STIFF WING 
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Mode 

Natural frequency 
(Hz) 

Reduced natural frequency 

1 

0.44 

0.0072 

2 

1 .20 

0.0198 

3 

2.33 

0.0385 

4 

4.03 

0.0666 



FIGURE 6 NATURAL MODES AND FREQUENCIES FOR THE FLEXIBLE WING 
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3.2 Forced Vibration: Frequency Response Function 


The final form of the equations for the frequency response 
function of the wing tip deflection due to a sinusoidal gust have been 
obtained by separating variables, integrating the spanwise influence, 
and assuming sinusoidal form for the solution. Remaining yet unsolved 
is a system of linear equations for the amplitudes of the different 
modal responses. The system of equations is: 




1 ^ r 2 1 2bpa(y*) 

J P„C„ - E [k\j - k C(k) B„jJ Cj = A K(k) 

w ' 


o) • b_ 

where = — y — ; reduced natural frequency of mode 


(46) 


fti(|),-'|)iCly/TrpbnS ; normalized general mass 
11 K 


-a 


A. . = 

ij s 


2 

a (J)^.(})jdy and B^.^. = 


-I 


ac})^.(j).dy ; normalized 


aerodynamic cross terms 


The unknowns, 5. = n^/boCV/u), are normalized amplitudes of the 
modal wing tip deflection in terms of the reference semi -chord, b^, and 
the magnitude of the disturbing force, a = IT/V. The remaining terms are 
the reduced forcing frequency, n., and the gust location, y*, which for 

vJ 

this study are spaced evenly at 38 locations along the span. The above 
system of equations is solved by Gaussian elimination for each forcing 
frequency and gust location. Fortunately, the coefficient matrix does 


25 


not change for each gust location, and the coefficient matrix is reduced 
only once in subroutine GAUSS for each forcing frequency. The nonhomo- 
geneous part is reduced for each gust location and the amplitudes are 
determined by back substitution in subroutine BACKS. 

Careful examination of the cross terms. A., and B.., shows that 

I J I J 

they vanish when the integration is of the product of an even mode and 
an odd mode. This means that the linear system of equations becomes 
uncoupled between even and odd modes, and the response of the system to 
a sinusoidal input as separated into response of even modes and response 
of odd modes can then be examined. Obviously, the response of the even 
modes to a gust located at y is the same as the response to a gust 
located at -y. The response of the odd modes to a gust is antisymmetric. 
In other words, the response of the odd modes to a gust at y is the 
negative of the response at -y. The system of linear equations is 
solved for gust locations on only half of the wing. Shown in Figures 
7-11 are the plots of the frequency response function versus the 
reduced frequency of the gust for several gust locations along the 
nonuniform wings. 

3.3 Spectrum Equation: Wing Tip Velocity Power Spectrum 

The final form of the output power spectrum due to a continuum 
of stationary, homogeneous, isotropic input is obtained from Equation 
(11), which can be integrated numerically using the trapezoidal rule to 
yield: 
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Reduced frequency -y— 


FIGURE 7 


Gust location at y (0.86 ft) 


FREQUENCY RESPONSE FUNCTIONS DUE TO GUST EXCITATION AT 
0.86 FEET FROM MIDSPAN 



Tip velocity of wing normalized by gust velocity 


FIGURE 8 


CM 



FREQUENCY RESPONSE FUNCTIONS DUE TO GUST EXCITATION AT 
14.76 FEET FROM MIDSPAN 
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Tip velocity of wing normalized by gust velocity 


O 


FIGURE 9 



Gust location at y (16.5 ft) 


FREQUENCY RESPONSE FUNCTIONS DUE TO GUST EXCITATION AT 
16.5 FEET FROM MIDSPAN 
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Tip velocity of wing normalized by gust velocity 


<M 

cn 



FIGURE 10 FREQUENCY RESPONSE FUNCTIOMS DUE TO GUST EXCITATION AT 
23.44 FEET FROM MIDSPAN 


30 








where N is the number of gust stations and A is the gust station width. 
The spectrum program determines the wing tip velocity spectrum, there- 
fore, Z. must represent the velocity at the right wing tip due to a 
sinusoidal gust at station i. The frequency response function is 
defined in the program as: 

. NM 

Z = ^ E n. (48) 

“ j=l ^ 

where NM is the number of elastic modes considered. The summation 
includes the response of only the elastic modes; the response of the 
rigid body mode is not included in the summation because the navigation 
system, located at the airplane's center of gravity, should be able to 
subtract the motion of the center of gravity from the turbulence data 
taken at the wing tip. There is a difference between the motion of the 

center of gravity and the motion of the rigid body mode; the center 

of gravity motion includes the motion of the even modes at the center 
of gravity. An attempt should be made to filter out the elastic mode 
motion from the center of gravity motion before correcting the turbulence 
data because the elastic mode motion at the center of gravity can be 
180 degrees out of phase with the elastic mode motion at the wing tip 
and can therefore introduce a larger error in the turbulence data. 

This study will assume that the motion of the center of gravity measured 
by the aircraft navigation system has been filtered so that it contains 
only the rigid body motion of the airplane before correcting the turbu- 
lence data. In this study wing vibration is defined as the motion of 
only the elastic modes of the wing. 
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Note that the frequency response function is normalized by the 

gust magnitude, u, and therefore represents the velocity at the wing 

tip due to a unit sinusoidal gust. The input spectrum is taken as 

von Karman's cross spectrum function of atmospheric turbulence [3] and 

is normalized in the program by the root mean square of the turbulence, 
2 

. Therefore, the wing tip velocity power spectrum is normalized by 
the root mean square of the turbulence. Reference [2] gives further 
examples of spectrum analysis problem solving. 
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4.0 RESULTS AND CONCLUSIONS 

A study of the frequency response function will aid in under- 
standing the power spectrum. An examination of the frequency response 
in the high frequency (reduced frequencies above one) domain displayed 
in Figures 7-11 shows that regardless of gust location, the wing tip 
velocity vanishes rapidly with increasing gust frequency. Physically, 
the wing is expected to respond less to higher frequencies because the 
higher modes have low response amplitude. Another gust location 
independent trend in the frequency response function is the extreme 
maximum at low frequency, also displayed in Figures 7-11. This is 
physically explainable because the gust frequency becomes close to the 
natural frequency of the wing. 

The more interesting trends of the frequency response functions 
of the wing are the gust location dependent trends. Figures 7-11 show 
the frequency resonse functions for both the standard and the flexible 
wing excited at five different locations along the wing. The locations 
of the gust are marked along the abscissa in Figures 5 and 6. Figure 7 
shows the response of the wing due to a gust near the midspan. 

Because the odd modes have their nodal point at the midspan, they 
should not respond as much as the even modes. The response is particu- 
larly clear in the curve for the standard wing where only the first and 
third elastic modes have large peaks. The response of the flexible wing 
is not quite as clear as the response of the standard wing because its 
natural frequencies are so close together the curve tends to be washed 
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out. Figure 8 shows the response of the wing due to gust excitation at 
14.76 feet from midspan. This gust location is at the maximum of 
second elastic mode deflection and close to the nodal points of the 
first and third elastic modes. The response curve should be like that 
shown for the standard wing. The frequency response curve for the 
flexible wing shows a much more pointed maximum than the response curve 
for the flexible wing shown in Figure 7 because in Figure 8 only the 
second elastic mode is participating while in Figure 7 the maximum 
consists of both the first and third elastic mode responses. The first 
elastic mode has its nodal point approximately 16.5 feet from the mid 
span. Figure 9 shows the frequency response function due to gust 
excitation at this point. The response of the flexible wing shows a 
pointed maximum characteristic of low first elastic mode participation, 
while the third mode shows a great deal of participation since the 
excitation is at its maximum deflection for the mode. As expected, the 
standard wing shows a very small peak for its first mode response. 
Figure 10 shows the response curve due to gust excitation near the 
nodal point for the second elastic mode. The response of the flexible 
wing shows a local minimum for the second elastic mode response, while 
the first and third elastic modes show peaks for their response. The 
standard wing shows no response around the second elastic mode. 

Figure 11 shows the response curve for gust excitation at 26.92 feet 
from mid span, which is near the nodal point for the third and fourth 
elastic modes. The curve shows minimum participation of the third and 
fourth modes for both standard and flexible wings. 
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The wind velocity measured at the wing tip is the sum of the 
turbulence, Up(t), at the wing tip and the velocity of the wing tip, 
U^(t), due to wing vibration. The measured wind velocity is then: 

U^(t) = Up(t) + U^(t) (49) 


The correlation function of the measured wind is: 

r T 


R (t) = -L 
T-^ 2T 


U„(t) V„(t + t) dx 
m m 


-T 


1 im 1 
T^ 2T 


T 

[Uw(t) U^(t + ,t) + Up(t) Up(t + t) 


J-T 

+ U^(t) Up(t + t) + Up(t) U^(t + t)] dx 

= Rw(t) + Rp(t) + R^p(t) + Rp^(t) 


(50) 


The Fourier transform of the correlation equation gives the 
spectrum equation 


d) = cb + d> + d) +,d) 

^m ^w ^p ^wp ^pw 


Since 4>^p = (J)p^, the equation can be simplified: 


<f) = (}) + (j) + 2Re(j) 

^m ^w ^p ^wp 


(51) 


(52) 


where c|)p is the power spectrum of turbulence at the wing tip and is 
the power spectrum of the wing tip velocity. Figures 12-15 show 
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Normalized power spectrum ((j> and (j> ) 

p w 



(ob 

Reduced frequency -y— 


FIGURE 12 SPECTRUM OF ATMOSPHERIC TURBULENCE FOR LENGTH SCALE 

132 FEET AND SPECTRUM OF WING TIP VELOCITY FOR FLEXIBLE 
AND STANDARD WINGS 
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Normalized power spectrum (t()p and (\,^) 



FIGURE 13 SPECTRUM OF ATMOSPHERIC TURBULENCE FOR. LENGTH SCALE 

2112 FEET AND SPECTRUM OF WING TIP VELOCITY FOR FLEXIBLE 
AND STANDARD WINGS 
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Normalized power spectrum (4> and (j) ) 


-Atmospheric turbulence 



FIGURE 14 SPECTRUM OF ATMOSPHERIC TURBULENCE FOR LENGTH SCALE 
132 FEET AND SPECTRUM OF WING TIP VELOCITY FOR STIFF 
WING 
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Normalized power spectrum (<t) and cj)^) 


FIGURE 



15 SPECTRUM OF ATMOSPHERIC TURBULENCE FOR LENGTH SCALE 
2112 FEET AND SPECTRUM OF WING TIP VELOCITY FOR STIFF 

W 


the wing tip velocity power spectrum and, for comparison, the atmo- 


spheric turbulence power spectrum. The cross spectrum (jj^p contains 

phase information between the wing tip velocity and the atmospheric 

turbulence velocity. The equation for 4 used in the program is: 

wp 

f 

<l>wp = Ua - S) ds (53) 

■'0 


The relative error introduced due to wing tip velocity is: 



(54) 


Figures 16-21 show this relative error as a function of frequency 
for the three different wings and different turbulence length scales. 

The flexible wing shows a large relative error throughout all the 
lower frequencies. An airplane having these wing characteristics is 
not recommended for measuring atmospheric cross spectra across its wing 
span. The standard wing shows a large relative error close to its first 
natural bending frequency. The first natural frequency can be close to 
the lower frequencies of atmospheric turbulence (approximately 0.04 Hz) 
if the turbulence length scale is small enough. For large turbulence 
length scale, the maximum error introduced due to wing vibration of the 
standard wing is reduced but still amounts to some 50 percent more over the 
turbulence spectra scale as uL/V and; hence, for higher length scales the 
frequency of the wing vibration is out of the frequency range which con- 
tains significant turbulence energy. If the wing has characteristics 
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Relative error 







Relative error 



o)b 

Reduced frequency 


URE 17 RELATIVE ERROR IN MEASURING TURBU 
,2112 FEET FOR FLEXIBLE/WING 



Relative error 









Relative error 







of the standard wing. It is recommended that accelerometers be mounted 
on the wing tip to obtain accurate measurements of atmospheric turbulence. 
The stiff wing shows a smaller relative error than the standard wing and 
the range of larger error is slightly shifted to the high frequencies. 

The relative error still reaches 50 percent and for small turbulence 
length scale can get close to the lower frequencies of atmospheric 
turbulence. If accurate measurements of atmospheric turbulence fre- 
quencies on the order of 0.22 Hz are desired when using an airplane with 
the stiff wing characteristics, it is recommended that accelerometers be 
mounted on the wing tip. This study concludes that to measure atmospheric 
cross spectra across an airplane wing span, a stiff wing is required and 
that to measure accurately the whole range of the atmospheric turbulence 
spectrum, accelerometers mounted on the wing tips are required. 
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APPENDIX A 


USER'S GUIDE 

This appendix contains details of the computer code used for 
this study. The four flow charts illustrate the complete computer code, 
the free vibration subroutine, the forced vibration subroutine, and the 
spectrum analysis subroutine. Following the charts is an explanation of 
how to modify the computer code for different aircraft. Finally, the 
complete computer code listing is provided. 
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A.l FLOW CHART OF THE COMPLETE PROGRAM 


NC-Number of Mode To Start On 
N ~ Total Number of Mode To Determine 
NN ~ Number Modes To Determine Mode On 
W ~ First Guess of Frequency 


I 


Free Vibration Program 


I 


W^ ~ Frequency 
M^- ~ Generalized Mass 
4)^- ~ Mode 

A^-j & ~ Aerodynamic Cross Product 


I 


Forced Program 


I 


k ~ Driving Frequency 
Z-j(y) ~ Response Of i Mode To Gust 


3 Input/Output 
H Program 


I 


Spectrum Program 


1 


k ~ Frequency 

(t)^ ~ Power Spectrum Of Wing Tip Velocity 
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A. 2 FLOW CHART FOR THE FREE VIBRATION PROGRAM 



I I Subroutine 


|yv 


Variables Passed 
Between Subroutines 
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A. 3 FLOW CHART OF THE FORCED VIBRATION PROGRAM 



(2) Input/Output 
I Subroutine 

I Variables Passed Between Subroutines 
IRK 
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A. 4 FLOW CHART OF THE SPECTRUM PROGRAM 



Q Input/Output 


I Subroutine 

|rK Variables Passed Between Subroutines 
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A. 5 PROGRAM MODIFICATION 


The computer code is actually three programs: the free vibra- 

tion program, run on the IBM 360 at The University of Tennessee, 

Knoxville; and the forced vibration and spectrum analysis programs, 
both run on the PDP 11 at The University of Tennessee Space Institute. 

The input-output data format is compatible between programs, and all 
three programs may be changed and extended. Listed below are changes 
that must be made to each program for different wings. 

The free vibration program was written for application to any 
wing. To adapt this program to a different wing, the function El, 
specifying bending stiffness along the span, must be changed. Also, the 
functions EIP and EIPP, which are the first and second derivatives of 
El, might need to be changed if El is more complicated than a step func- 
tion. The array RA, which defines the wing's semi-chord along the span, 
and the array X, which defines the mode position, must be changed. 

They are defined in subroutine STRK. Finally, the function RM, the mass 
per unit length, must be corrected. 

The forced vibration program is extremely simple to modify for 
other aircraft. The array X, specifying the gust locations, must be 
changed in MAIN. Also, in MAIN the constants BR, reference semi -chord; 

S, surface area of wing; and U, mean flight speed; must be changed to 
fit the airplane. Finally, the function RA, wing semi -chord distribution, 
must be altered. 

There are only two constants that must be changed to modify the 
spectrum analysis program for other airplanes. These are the reference 


55 



semi -chord, BR, and the mean flight speed, U. This program is most 
likely to require change due to a new distribution of the atmospheric 
turbulence cross spectrum, which will require that function TSPEC be 
rewritten. The program can be modified to give the power spectrum of 
the root bending or wing deflection at any point of the span by 
redefining the response function, the array Z. 
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A. 6 PROGRAM LISTING 


CaMMON/EOAT/EEE 

EEEsBOOQQOOOOO^O 

CALL D1 

CALL 02 

CALL 03 

STOP 

£l^D 


SUnKUUTINE 01 

C*** THIS PfiUGRAM IS USEO FOR OCTERMlNti WING'S NATURAL 
C FREOUENCI AND MODE, 

C 

C IT USES A BISECTION LIKE METHOD TO IMPROVE THE GUESS 
C ON THE FREtiUENCy AND MAKES THE CHARACTERISTIC DETERMINANT 
C VANISH. 

C 

C RUNGE-KUTTA FEHLBERG 7L8 IS USED TO DETERMINE SOLUTION 
C TO DIFFERENTIAL EON IT JIAS VAIRABLE SIZING BETWEEN 
C FIXED NUDES. 

C 

C SlilPSOUS METHOD IS USED TO INTEGRATE FOR THE AREOOYNAMIC 
C CROSS TERMS AND GENERALIZED MASS. 

C 

C THE PROGRAM ilAYbE ADOPTED BY CHANGING FUN...SEI 

C (TilE FUliCTlON DESCRIBING THE DISTRIBUTION OF THE PRODUCT 

C YOUNG'S ELASTICITY AND SECTKjnAL mASS MUENT, IF 5E1 IS TO 

C HAVE DERIVATIVES UP TU THE SECOND ORDER FUN...FUNEV NEEDS TO BE 

C CHANGED TO INCLUDE THEM IN THE DIFFERENTIAL EON. THE WING 
C plan ARRAY, RA, WILL NEED T0 BE CHANGED, 

C** + NC = TllE number of mode TO START CALCULATING 
C**+W=THE FRIST GUESS OF THE NATURAL FREQUENCY OF THE NC MODE 
C**+H=STEP SIZE FOR SEARCH ROUTINE 
C***N=TutAL NUMBER OF MODES TO CALCULATE 
C***NN = iJuMBER OF NUDES TO CALCULATE MODE ON 
REAL*8 W,H 
NC = 2 
H=2,0D0 
W=4.0D0 
N = 4 

Ni'J = 151 

CALL STRKCNN) 

CALL SZERO(W,H,N,NC,NN) 

RETURN 

END 
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SUBROUTINE S^ERO ( W , H , N , NC , NU } 

C THIS subroutine STEPS AUUNG THE FHEQUENOf DOMAIN ONTIDL 

C IT bRACKEfS THE NATURAL FKLOUENCY AND THE CALLS A BISECTION 

C ROUTINE AND A SECANT ROUTINE Tu IMPROVE THE WIDTH OK THE BRACKET 

C The SUbRUUTINE FINALLY CALLS SUB...DMUDE TO DETERMINE THE 

C MODE, AND THEN CALLS COEF TO SET UP THE INTERGRATIONS OF MODES 

c And their products. 

C**** EL E EUsLRROR BOUNDS FOR KUNGE KUTTA 

C**** ril L W2=BKACKETS FOR FREQUENCY AND STEPS FOR THE SEARCH ROUTINE 
C**** FI & F2=VALUES OF THE CHARACTERISTIC DETEKMINTE FOR W1 S. W2 
C**** N=NUMBEK Of EIGENVALUES TO BE SEARCHED 
C**** W=GUESS FOR EIGEN VALUE 
C**** H=oTEP SIZE FOR SEARCH 

C**** i\,C=NUMBEK OF THE FRIST MODE TU CALCULATE 
C**** NN=NUMBER OF NODES TO CALCULATE MODES ON 
C**** l=NUMbER OF MODES CALCULATED DURING PROCESS 
REAL*8 W.H.Wl ,W2,F1 ,F2,EL,EU,Y,X 
CUWMUN/FACT/Y(8,15U ,X(151) 

DIMENSION YYU51) 

1 = 0 
W1 = W 

EL=.0001D0 

EU=.001D0 

103 CALL FUNCWl ,F1 ,EL,EU,NC,NN) 

101 W2=ttl+H 

CALL FUN(W2,F2,EL,EU,NC,NN) 

1FIF1+F2.LT.0) GO TU 102 
W1 = W2 
F1=F2 
GO TO 101 

102 CONTINUE 

CALL B1SEC(W1,F1,W2,F2,NC,NN) 

CALL SECAN'KWl ,F1,W2,F2,NC,NN) 

CALL DMUDE(YY,NN) 

WRirE( 6 , 200 ) W2 
200 FUK.MAT(2X,F20. 10) 

WKil E{b, 200 H y Y( J) ,J=1,NN) 

1 = 1 + 1 
«S=W2 

Call cof(Nc,yy,ws,nn,I) 

NC=NC+1 

IFCl.GE.N) GO TO 104 

W1=W2+H 

GO TO 103 

104 RETURN 
END 


SUBROUTINE FUN ( W , F , EL , EU , NC, NN) 

C THIS SUbKtPUTIIME FIXES THE INITIAL VALUE FUR THE SOLUTION 

C Al'lD THEM CALLS SUb . . . RK7 C RUNGE KUTTA ROUTINE) TO DETERMINE SOLUTION 

C AFTER WHICH THIS SUB CALCULATES THE VALUE OF THE CHARACTERISTIC 

C DETERMlNAiiT. 

C*** EL t, EU =ERROR BOUNDS hOR RUNGE KUTTA 
C*** ,'JC=WHlCtl MODE WORKING ON DETERMiNG 

C*** Y=FUNDAMENTAL SOLUTIOWS TO DIFFERENTIAL EQN OUTPUT FROM RK7 
C*** F=VALUE OF CHARACTERISTIC DETERMINTE 
REALtB RK,W,F,EL,EU,Y,X 
COMWUN/FACT/YC8, 151) ,XU51) 

COMMON/EGnVL/RK 
KK = W 
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n n n n n r. n 


C INITIAL COriOlTlON TEST FUR EVEN OR ODD MODE 
DO 100 1=1,8 

100 Y(1,1)=O.DO 

lF(NC.ECi.2.0H,NC.EQ.4} GO TO 101 - 
C INITIAL CONDITION FOR ODD MODE 
YC2,1)=1.D0 
YC8,1)=1.U0 
GO TO &00 

C INITIAL CONDITION FOR EVEN MODE 

101 YU, 11 = 1, DO 

YC7,1)=1.D0 

500 CALL RK7(8,NN,EL,EU) 

F=YO,NN)*Y(8,NN)-Y(7,NN)*Yt4,NN) 

RETURN 

END 


SUHKiJUTIWK BISEC(X1 ,F1,X2,F2,NC,NN) 

REAL*8 XI ,X2,F1 ,F2,FH,XH,EL,EU 

CCB=.0l 

EU=.001D0 

EL=.0001D0 

102 X.M=(X1 + X2)/2,D0 

CALL FUNCXH,FM,EL,EU,NC,NN) 

lF(Fl*fM.LE,0.D0) GO TO 100 

X1=XM 

F1=FM 

Gu TO 101 

100 X2=XM 
F 2=Ff1 

101 RE=DABS(X1-X2)/X1 

IF (KE.GT.CCB) GO TO 102 

RETURN 

End 


SUBROUTINE FUNEV C K , X , Y , F) 

THIS SUllRUUTlUE SUPPORTS THE RUNGE KUTTA AND DESCRIBES THE 
DIFFERENTIAL EON 

THIS SUB WILL NEED CHANGING IF 5EI IS TO HAVE DERIVATIVES 
OR IF THE T<-1RS10NAL MODES ARE BEING DETERMINED 
♦ K=ORDEK ()F THE TAYLOR SERIS TERMS 
**** F=DER1VAT1VLS VALUES 
C**** Y=F UwDAHEIJTAL SOLUTION 
C**** W=FREOUEWCY 

REAL+8 F, Y ,DRM,W,SE1,X 
DIMEWSION FC8,13),Y(8) 

CIJMi'inil/EGNVL/W 
f U,K)=Y(2) 

F(2,KJ=Y(3) 

FC3,M=Yt'H 

F(4,M=DRM(X)*W**2*Y(1)/SEICX) 

F(5,K)=Y(bJ 
F C 6 , K ) =Y ( 7 ) 

F(7,K)=Y(8) 

F(8,M=DkM(X)*W**2*YC5l/SEI(X) 

REiURIJ 

END 
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FUNCTIUN SEUX) 

C*** X=DiaTAMCE KKOrt SEMI SPAN 
l/'iPLIClT RKAL + 0 (A-H,C-ZJ 
SE1=3000000000,DO 
IF(X.GC.n.UO) SE1S300000000.UO 
RETURN 
END 

SlJbRUUTlNE SECAMTCXl ,F1,X2,F2,NC,NN) 

Kt.A)j*8 XI ,F1 ,X2,F2,XHl,t'P,OX,XM,-fH^CCS,EU,EL 

Ei,= .ooooooooiuo 

EU=. 000000011)0 

CCb= *000)1)0 

XH1=X1 

CALL KUNCXl ,F1 ,Eli,EU,NC,NfJ) 

CALL FUi'l(X2,F2,EL,EU,NC,NN) 

103 FI-=CF2"Fl)/CX2-Xn 
UX = -t'l/l'P 
XH=X1+DX 

CALL FUl'j(XM,Fr'l,EL,EU,NC,hN) 

WRiXt;C6,200JXl ,F1 ,X2,F2 

200 FURMATC2X,2(D21.1A,2X,D21.14,2X)) 

C IF(DABS(XH-XM1)/XM).LT.CCS) GO TO 101 
1F(DAB5(FM).LT.CCSJ go to 101 
It (FM+Fl.LE.O.DO) GO TO 500 
X1=XH 

FI = fM 

XM1=XM 
Go TO 103 
500 X2=XM 
F2=t M 
X.‘11=XM 
GO TO 103 
101 X2=XM 
F2=FM 

'ArtITb(6,201) X2,F2 

201 FtJRMAT(2X,D21.14,2X,D21.14) 

RETURN 

END 

SuaROUTlNE UMODECYY ,UN) 

KlAL+H Y,X,C1,C2,D 
CUMMUM/FACT/Y(8,151) ,XC151) 

DlrtENSJiJN YY(15U 

D=Y(1 ,NN}*Y(7,NN)-Y(5,NN3*Y(3,NN) 

ClsYC7,NN)/r) 

C2=-YC3,NN)/D 
DU 100 1=1, NN 

100 YY(1)=C1*Y(1,1)+C2*Y(5,I) 

RETURN 

END 

function ORM(X) 

Implicit heal*8 (a-h,o-z) 

DRM=130.DO 

lF(X.Lt:.4.D0J DRM = 2205,D0 

IFCX.LE. 1 1 .D0.AND.X.GE.8.D0) DRM=2600.D0 

RETURN 

End 

function RMCXJ 
KM=130. 

IFCX.LE. 4.) RM=2205. 

IF(X.LE.11..AHD.X.GE.8.) RK=2600. 

RETURN 

END 
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FUNCTION S1HP(Y,H,N) 

DIMENSION Y(NJ 
T1=0. 

Jl=N-2 

T2=0. 

DU 100 1=3, Jl, 2 

100 Tl=Ti+Y(l) 

DO 200 1 = 2, N, 2 
200 T2=Y(1)+T2 

SIMP=H*(y(l)+Y(N)+2.*Tl+4,*T2)/6, 

RETOKh 

END 

SUBROUTINE COF(N,YY,w2,NN,Il) 

C THIS SUBROUTINE SETS UP THE FUNCTIONS FOR INTEGRATION 

C*** YY=NEW MUUC THEN LATER USED AS SCRATCHED ARRAY 
KP=AKKAY OF MODES 
C*** NsFREOUENCY 
C*** Gi'1=GENEl<ALIZED MASS 
C*^* A «. b = AKEi)DYNAMMIC CROSS TERMS 
c**^ ra=vving plan 

C**» fJsNUHBER OF MODE BEING WROKED ON 
C*** NW=UUMhER OF NUDES 
C*** 11 =imUMBER of TIMES COF CALLED 
CUMHON/DATARA/RA (151) 

COHMuN/TRAN 1/W(5) ,GMC5),RP(5,151),A(5.5),BC5,5) 
UlMEiiSIUN YY(15U,XX(151),RR(151) 

C*** READ DATA 

H=3i,/(tJti-l ) 

1F(I1.<;t. 1) GO TO 500 
DO 601 1=1, NN 
bOl XXCI)=(1-1 )*H 
uM = iJ-l 
Wl 1 J=0, 

GMC 1 )=40000.0 
DU 600 1=1, NN 
600 HP(1,1)=1, 

AC 1 , 1 )='15. 036925 
iUl,l)=52. 939763 
C*** CALCULATE NEN DATA 
500 CU1<T1i4UE 

DU lOl 1=1, NN 

101 RP(N,I)=yY(I) 

DO 102 1=1, NN 
A1=K.'UXX(1) ) 

B1=RP(N.1)*RP(N,I) 
liYU) = Al<=Bl 

102 CONTINUE 

G/1 (N)=SIMP(YY,H,NN)*2. 
w(N)=W2 
DO 104 d=l,N 
NA2=N/2. 

NRl=J/2. 

H2=N-7*HR2 

Nl=J-2.«uRl 

IF(N2.EQ.0.AND,ia,EO,0> GO TO 400 
IFtJ.EO.l) GO TO 400 
lF(N2.E<J.l.AND.Nl,Ea,l) GO TO 400 
A(N,d)=0.0 
B(N,J)=0.0 
GO TO 105 
400 CUNIINUE 
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n f 


103 


DO 103 1=1, hN 

T(Ul)=RAtI)*RPCd,l)*RPCN»I) 

RR(I)=RA(l)*yY(I) 

A(N,J)=SIMP(RR,H,NN) 

B(N,0)=&lMP(yy,H,NN) 

105 CUNllriUE 
104 CONTINUE 
RETURN 
END 

SUBROUTINE STRKCNNJ 

THIS SUBROUTINE INITIALIZES THE COEFICIENTS FOR RK7 
FIXES THE NOUES AND WING PLAN 
C**** Ni4 = NUMtiER OF NODES 
C**** A,B,C = C0EFFICINX6 FDR RK7 
C**** X=NUDE VALUES 
C**** HA=WING PLAN ARRAY 

REAL*8 A,B,C,CH,X,y ,DEL 
CUMMON/UATARA/RA( 151) 

CUHMUN/RKC/AU3) ,BC13,12),C(13),CH(13) 
COMMUN/FACT/y(8,151),X(151) 

DEL=33.D0/(NN-1) 

DEL2=33./(NN-1D 

ND=CNN-l)/3, 

RMl = (tjU-U /a. 

Du 100 1=1, NN 
X(I)=tI-l )*DEL 
KACI)=1.0 

IFU .(;r.ND)RA(I)=l,-,027*(I-RMl)*DEL2 

100 CUNl'ILllE 

ACl J=0.0D0 

A(2)=2. DO/27. DO 

A( J)=l .DO/9.DO 

A(4)=1.D0/6.U0 

A(5)=b.D0/l 2.D0 

Alfe)=l .00/2.00 

AC7)=5.O0/6.D0 

A(8)=l .O0/6.D0 

AC9)=2.D0/3.D0 

AC10)=l .D0/3.D0 

AC11)=1.D0 

AU2)=0.D0 

AC13)=1,D0 

Bll.l)=O.0D0 

B(2, I )=2.D0/27.U0 

BC3, 1)=1. DO/36. DO 

B(3,2)=l. DO/12. DO 

BC4,1)=1.U0/24,D0 

B(4,2)=O.DO 

B(4, 3)=1 .D0/8.D0 

B(5,l)=5. DO/12. DO 

BC5,2)=O.ODO 

BC5.3)=-25. DO/16. DO 

Bl5,4)=25. DO/16. DO 

B(b, 1 )=1 .D0/20.DO 

b(6,2)=0.D0 

BC6, 3)=0.D0 

B(6,4)=1 .D0/4.D0 

B(b, 5)=1 .D0/5.D0 

BC7, 1 )=-25.D0/108.D0 

D(7,2J=0.00 

B(7, 3)=0.0D0 

Bt7,4)=125.D0/108.D0 

B(7,5)=-fa5.D0/27.D0 
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B(7,6)=12b, DO/54, DO 
B(8, n=31 .DO/30O,DO 
bCti,2)=0.0D0 
BC8,3)=0.D0 
B(8,4)=0.D0 
B(8,5)=bl.DO/22b.OO 
B(8,6)=-2.D0/9.U0 
BC8,7) = i3, 1)0/900, DO 
BC9, n=2.D0 

B(9,2J=O.DO 

6(9, 3)=0.D0 

B(9,4)=-53,00/6,DO 

B(9,5)=704.D0/45.D0 

BC9,6)=-107. DO/9, DO 

B(9,7)=67.U0/90,D0 

D(9, 0)=3.D0 

B(10,1)=-91,D0/108.D0 

B(10,2)=O.ODO 

B( 10 , 3 7=0 .DO 

b( 10,4)=23.D0/108.D0 

B(10,b7=-976,D0/135,U0 

BCIO, 07=311, DO/54, DO 

B(lO,7)=-19,D0/60,DO 

BCIO, 87=17, DO/b, DO 

B(10,97=-l. DO/12. DO 

BUI , 17=2 3«3,D0/41O0,D0 

b( 1 1 , 27=O.DO 

b( 1 1 , 37=O.LiO 

b(ll,47=-J41. DO/164. DO 

BUI , 67 = 4 496.00/102 5.00 

B(11,67=-301. DO/82. DO 

BUI, 7 7 = 2 1 33, DO/4 100, DO 

BUI , 87=45. DO/82. DO 

HUl, 97=45. DO/164. DO 

BUI, 107 = 18. DO/41. DO 

b(12,17=3.D0/205.D0 

b(12,27=0,D0 

B (12,37=0. DO 

b(12,4)=O.DO 

b (12, 5 7 =0.1)0 

BU2,67=-6. DO/41. DO 

BU2,77=-3,D0/205.D0 

13(12, 87=-3. DO/41. DO 

BU2, 9 7=3.00/41 .DO 

b(12, 107=6.00/41.00 

BU2, 117 = 0.000 

B(13,17=-1777.D0/4100.D0 

BU3,27 = O.DO 

flU3.37=0.D0 

BU 3, 4 7 =-341 .DO/ 164. DO 

B(1 J,b7=4496.D0/io25.D0 

u (13, 6 7 =-289. 00/ 8 2. DO 

B( 13, 7 7 =2193. 00/4100. DO 

B(13, 87=51. DO/82. DO 

B( 13,97=33.00/164.00 

BU3, 107 = 12. DO/41. DO 

0(13,117=0.00 

BU3, 127 = 1.00 

CU 7=41 .D0/840.D0 

C(27=O.DO 

C(3)=0.00 

C(47=O.UO 

C(57=0.00 
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C(6)=34.D0/105.D0 

Cl7)=9.O0/3b.D0 

CC8J=9. 1)0/35, DO 

C(9)=9.u0/280,D0 

C(10)=9.U0/280.D0 

CCin=41.D0/840.D0 

C(12)=0.l)0 

C( 13)=0.U0 

CHll)=0.00 

CHC2J=0.1)0 

C)U3J=0.D0 

Cri(4)=O.UO 

CH(5J=0.U0 

CH(t>) = 34. 00/105, DO 

CH(7)=9, DO/35, DO 

CH{8)=9, DO/35, DO 

CHC9)=9.D0/2H0.D0 

Cfi( 10)=9,D0/2b0,n0 

CHCin=0,D0 

Crit 12)=41 ,D0/840,D0 

Cil(131=41,D0/840.D0 

RETURN 

SUBRUUTlNt; KK7(Nf.,NN,EL,EU) 

RUilGh; KUTTA FEHDBERG StVhNTH ORDER 
C*** el=c;rk(jh duwek dui.iwu 

C*** El) = hKKOK UPPER DlllIND 
C*** NS = l\iUMMER OF SYSTEM OF EON 
C*** Y=.-itJDUl’10N 

N0 = IJUMBER OF PTS TO DETERMINE THE SOLUTION 
C*** RL=LENGfH OF INTERVAL 

DIMENSION Yo(8) ,F(8 ,13) ,YY(8) 

DIMENSION DY4(0),DY5C8),Y1(8) 

CUMMOn/KKC/A( 13) ,U(13,12),Ca3) ,CH(13) 

CumM0N/FACT/Y(8, 151) ,X(151) 

REAL* 8 Dy4,DY5,XX,YY,TH,A,B,C,F,DDl ,DD2,H,Y1,Y0,CH,EL,EU,X,Y 
DO 101 1=1, NS 
101 Y0U)=YC1,1) 

i><T = NN-l 
L=1 

C*** MAIN DU LOOP INCREAMENT TO EACH NODE 
DU 100 11=1, NT 
NC=0 

H=x(ii+n-x(ii) 

Gli to 2 03 

207 L=L-1 

Go TO 203 
206 L=L+1 

NC = 1 

203 DU 201 1=1, NS 
201 Yl(I)=YO(I) 

TH=H/L 

C*** DU LOOP FOR STEPS BETWEEN NODES 
DO 200 12=1 , L 

C*** DETERMINE THE NEEDED FUNCTION EVALURATION 
DU 300 K=l,13 
Kh=K-l 

DU 301 J=1,NS 

XX=X(I1 )+TH*C12-l)+A(K)*TH 
YYtJ)=Yl(J) 

IF(KM.EO.O) GO TO 303 
DO 302 13=1, KM 

302 YYC<J)=TH*0(K,I3)*F(J,I3)+YY{J) 
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303 CUIJTINUE 
301 COW'flNUfc, 

300 CALI. FUNEVCK,XX,YY,F) 

C*** UETt-KMItiE SOLUTION VALUE FOR END OF STEP 
DO 500 1=1, NS 
UY4(I)=0.0 
500 DY5(I)=0.0 

DO 401 1=1, NS 
DO 402 K=l,13 

DY4(1)=1'H*C(K)*F(I,K)+DY4(I) 

402 DY5C1 J=TH*CH(K)>t!F(I ,K)+UY5(I) 

401 CONTINUE 

c*** error and step size control 

DDl=UAbS((UY4(U-DY5(l))/(DY4(l)+Yl(l))) 

UU2=UADS( (DY4(5)-UY5(5) )/(DY4(5)+YlC5) )) 
1FCDU1.1 jT,£U.AND.DU2,LT.EU) go to 202 
GO 10 206 

202 1FCDD1,GT.EL,AND,DD2,GT,EL) go to 204 
IF(L.EO.l) GO TO 204 
IFCWC.EiJ.U GO TO 204 
GO TO 207 

204 CONTINUE 

DO 205 1=1, NS 

205 YUX) = yi(I)+DY4CIJ 
200 CONTINUE 

DO 102 1=1, NS 
YOCn = Yl (I) 

102 Y(1,I1 + 1)=Y(J(1) 

100 CONTINUE 
RETURN 
END 

SUBROUTINE D2 

C THIS PROGRAM DETERMINES THE AMPLITUDES OF THE Dlf FERENT 
C MODES TO A SINUSODAL GUST AT THE DIFFERENT STATIONS ALONG 
C THE WING. SUB... DO TAKES CARE OF INPUT AND OUTPUT PLUS 

C SETS UP THE CUEEFICIENTS THAT ARE DRIVING FREOUENCY INDEPENDENT 

C SUiiKUUTlNE COEf SET UP THE CUEFFICTEhT MATRIX FOR EACH DRIVING 

C FREOUENCY. WHILE SUb... GAUSS DOES HALF OF THE REDUCTION AND 

C SUB... HACKS FINISHES THE REDUCTION AND DUES BACK SUBSITUTION 

C FUR THE different NUN-HOMOGENOUS VICTORS CORRESPONDING TO 
C DlffEREtlT GUST LOCATIONS. 

C*** Nn=NUMHER UF NODES THE MOOES ARE DETERMINE ON 
C*** n2=NUMBER of gust LOCATION 
Nn=151 
N2 = 20 

CALL UO(NN,N21 

RETURN 

ENU 


SUBROUTINE 00(NN,N2D 

C THIS SUBROUTINE TAKES CARE UF INPUT AND OUTPUT 
C AND PERFORMS OPERATIONS THAT ARE DRIVING FREOUENCY 
C INuEPENDENT, MEANING GAMA AND OMEG, 

C*** Y=SOLUTION, AMPLITUDES OF MODES 
C*** RP=MODE array 
C + + * X=GUST LOCATKHI 

C*»» A,b=AERODYNAMlC CROSS PRODUCTS 
C + ** W=NAlUt<AL FREQUENCIES OF MO'UES 
C**4 GH = GEi'JERAL.17.ED MASS l)t MODES 
C*** OMEG=kEDUCED natural FREOUENCY 
C*** GAMA=NUN DIMENSIONAL GM 
CUMHON/SOL/Y(20,5J 

CUIiMUN/TRANl/W(5) ,GM(5) ,RP(5, 151} ,A(5,5) ,B(5,5) 
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CUMMUI)/TRAN2/KKK(37) 

CUK.MUU/TKAM3/SY(37,19,5) 

CUWftUN/tDAT/EEfc; 

CU(«lMUU/FAC/PI,bR,S,KO,U • 

DXMENSIiJN X(20) 

CUHPUtX Y,SY 
CUMPbEX CMPLiX 
CUMMUN/UAT/GAMA(5) ,0MEG(5) 

N=37 

C50 FURMAT(2X,E14.7) 

C Wk1'IE(7,S0) EEE 

Pl=3. 14159 
BR=19. 0/2.0 
5=960.0 
R0=.0765 
U=575.0 

C*** PEKFUKM ARITilMATIC 
DU 601 1=1,5 

GAMAd 1=GH(I)/(PI*H0*S*BR) 

UMfc;G(l)=W(l)*UR/U 
JN = I 

DU 602 J=1,JN 
A(l,J)=A(I,J)*bR/S 
bCl ,Jl=b(I,JJ*BR/S 
A(J,1)=A(I,J) 
b(J,l)=B(I,J) 

602 CONTINUE 
601 CONTINUE 

DO 500 1=1,37 
XF(l.LC.lO) RK=1/1000. 

IFCI.L.E.19.AHU.I.GT.10) RK=(l-y 1/100, 
1FCI.LE.28.AND,I.GT,19) RK=(I-18)/10, 

IF(1.GT.28) RK=CI-27) 

DEL=33,/(N2-1) 

N2HSN2-1 
DU 100 11=1, N2M 
100 X(U )=Dbli*(Il-,5) 

CAl.L CUEF(RK,X,DEb,N2M,NN) 

CALL GAUSSt5,N2M) 

CALL BACKS(N2K,5) 

RRK(11=RK 
DU 1 J=1,N2M 
DU 2 J2=l,5 
SY(l,vJ,J2)=Y(J, J2) 

2 CUNTliiUrJ 
1 CUNTINUE 
C W1UTE(7,101) RK 

ClOl FQRHAT(2X,E14,7) 

C DU 102 J=1,N2M 

C VkiaTE(7,103HY(J,Il),H=l,5) 

CK 3 fOKHAT(lX,3(2E13.6)) 

C102 CUNTINUE 
500 CUNTINUE 
RETURN 

SUURUUTINE COEF(RK,X,DEL,N2H,NN) 

C TH15 SUBROUTINE SETS UP COEFFICIENT MATRIX AND THE DIFFERENT 
C NUN HUhUGEHOUS VICTORS. 

C** C=C0EFF1CIENT MATRIX 

C** 0= ARRAY UF NUN HOMOGENOUS VICTOR 

C*** C=CUEFICIENT ARRAY 

C*** RK=REDUCED FREQUENCY 

c*** x=lucatiun of gust 

CUMMUN/LS/C15.5) ,DC20,5) 
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C0MMUi'</DAT/GAMAC 5) ,0MtGC5) 

CUMMari/i'KANi/W(5} ,GM(S) ,RPCS*t5n*A(5,5).BC5,5) 

CUMMUU/FAC/P1,BR,S,«0,U 
Dlhbi^SlUM X(20J ,51(5,5) 

CUHFLKX C,CI,l),CC,RKK 
CurPBKX CMPliX 
C*** READ DATA 

C1=CMPDX(0.,1.) 

DO 103 1=1,5 
DU 104 J=l,5 

C(l, J)=-HK**2*A(I, J)+2*C1*RK*CC(RK)*B(I, J) 

IF(l.NE.J) GU TO 901 

C( X , J)=C(I , J)+GAMA(I)*(OMEGCI)*»2-RK**2) 

901 5U1,J)=CABS(CC1,J)) 

104 CUI-iTIIJUE 
103 CUNTIi'fUE 

C DO 902 1=1,5 

DU 106 J=1,N2M 
DU 105 1=1,5 

D(J, I )=2*BR/S*RA(X(J) )*RPH(I,X(J) ,NN)*RKK(RK)*DEL 

105 CUNTINUE 

106 CUi'JTIi'JUE 
KETURM 
END 

FUNCTION RPtl(I,Y,N) 

C0MMUN/TRAN1/WC5) ,GM(5) ,RP(5,151) ,A(5,5) ,B(5,5) 

UEL=33,/(N-11 
Ni'j=Art6(Y) /DEL>1 

R3=Ab5CY)/DEL+l,-NN 

KPII=KP(1,NN)+RS*(RP(1,NN+1)-RP(I,NN) )/DEL 
IF(Y.DT.O.O) GO TO 500 
GU TU 600 
500 CONTINUE 

IF ll.E0.3.0R.l,EQ.5) RPH=-RPH 
600 RETURN 
END 

FUiJCTIUN RA(Y) 

RA=1-,027*(AH5(Y)-11.) 

IF(AbS(Y) .LE.ll) RA=1, 

RETURN 

END 

FUNCTION RYICX) 

Z=(X/3. )**2 

RY1=( ( (( (( .0027873*Z-.0400976)*Z+.3123951)*Z-1.3164827)*Z 
++2. 1682709) *Z+. 221 2091) *Z-,6366198+. 63661 98*X*AL0G(X/2,)*RJ1(X) 
+ )/X 
RETURN 
End 

COMPLEX FUNCTION CC(RK) 

C THIS FUNCTION CALCULATES THE THEODOBSEN FUNCTIONS 
COMPLEX Cl 
CUHFLEX CMPLX 
C1=CMPLX(0.,1.) 

PJ1=RJ1 CRK) 

POO = RJO(KM 
PYI=Rlf 1 (KK) 

PYO=RYO(KK) 

F=PJ1*(PJI+PY0)+PY1*(PY1-PY0) 

G=PK1 ♦PYO+PJl+PJO 

H=(PJl+PY0)»+2+(PYl-PJ0)**2 

CCsCF+Cl^Gl/H 

RETURN 

END 
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COh'i^litX FUMCTION RKK(RK) 

C FUijCi’iUfi DtTEKHlNES THE GUST FORCE FUNCTION 
CUKHLiEX CI.CC 
COMPLEX CMPLX 
Ct=CMPLX(0. ,1,) 

PJ l = HJl (RK) 

RKx=CC(KK)*(RJO(RK)-Cl^PJl) +CI»P01 

RETORU 

E(mO 

FUNCTION RJl(X) 

Z=(X/3.)**2 

PJ1=C ( ( ( ( ( .0000 11 09+Z-, 0003 1761) *Z+, 00443319 )»Z-, 03954289 )*Z 
++, 21093573 )»Z-. 56249985) #Z+, 5) »X 

RETURN 

ENO 

FUNCTION RJO(X) 

Z= tX/3.)**2 

R00=( etc ( .000 21 *Z-, 0039444) *Z+, 044447 9) *2-. 3 163866 )*Z+ 1,2656208 )♦ 
+Z-2.2499997)*Z+1.0 
RETURN 
END 

FUNCTION RYO(X) 

Z=CX/3. )**2 

KYO=t ( C (C-, 00024846*Z+.00427916)*Z-. 04261 214) »Z+. 253001 17) »Z 
+-, 74350384) *2+. 6055936 )*Zj-, 3674669 1+,6366198*RJO(X)*ALOGCX/2.) 
RETURN 
END 

s 1) ti R 0 u.T 1 rj p: g a u s s ( n , n 2 ) 

C THIS SUiiRuUTli'JE DUES GAUSSIAN ELIMINATION FOR ONLY THE COEFFICIENT 
C MATRIX. SCALED PARTIAL PIVOTING IS USED, 

CuMililN/LS/CCS, 5) ,0(20,5) 

CUMMUN/PIVUT/IPEN (5) 

UIMENSIUH S(5) 

COMPLEX C.l) 

C*** K = PlViJT INDEX 
C*** C=CUFF1C1ENT ARRAY 
C*** 0 = lNituHfJGi-'.NOUS VICTOR 
C*** N = fJU;iHER OF EON 
DO 103 1=1, N 
IPENd )=I 
S(I)=0. 

DO 104 J=l,N 

104 IF(CAUSCC(I,J)).GT,S(I)) S C I ) =CABS ( C C I , 0 ) ) 

103 CONTINUE 

i>lM = i)-l 

DO 100 KK=1,NM 

I5=KK+1 

IP=IPEN(KK) 

J=KK 

CM=CAt)5(C(IP,KK))/S(IP) 

DO iOb 1=IS,N 
lP=lPEiUI) 

T=CABS(C(IP,KK))/S(1P) 

IF(T.LE.CM) GO TO 105 
CH=T 
J = 1 

105 CuNl'lNUL 
IPK=lPEiKJ) 

1PEN( J)=IPEN(KK) 

1PEN(KK)=IPK 

DO 101 1I=1S,N 

l = lPLricll ) 

K=1PF.N(KK) 

C(I,KK)=C(1,KK)/C(K,KK) 
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DO 102 J=1S,N 

102 CCl»d)=C(I,J)-C(l,KK)*C(K,J) 

101 CuNTlNUK 

100 COKTlhUE 
RETURN 
END 

SUHRUUTllJE BACKS(N1,N2} 

C DUbS kEDUCTIUN ON THE NON HOMOGENOUS VICTOR AND THEN DOES 
C BACK SUliSITUTIUN, 

N1=NUHHER UK NUN HUMOGENUUS VICTOR 
C*** N2 = L>lhE('iSlUN OK MON HOMOGENOUS VICTOR 
CUM MOM / K I VOIV 1 PEN ( 5 ) 

Coi'-.MUti/l,S/C(5, 5) ,DC20,5) 

CUMrtUi'J/S01j/Y(20,5) 

CUMPUEX c,D,y 
C**+ K1=SijUUT1(IN INDEX 

C**+ REDUCTION ON NUN HOMUGENOUS VICTOR 
DO 100 Kl=l,Nl 
1P=U>EN ( 1 ) 

Y(K1,U=1)(K1,IP) 

UU 101 KK=2,N2 
K = 1PEN CKK) 

•f = 0.0 

Jw=hK-l 

DU 102 J=1 , JN 

102 T=CIK,J)»Y1K1,J)+T 

101 YlKl ,KK)=DCK1 ,K)-T 

Y(Kl ,N2)=Y(K1 ,N2)/CCK,N2) 

C*** BACK .SUUSITUTION 
JJ = N2 

UU 103 K=2,N2 
US— UvJ 
JJ=JJ-1 
KKslPEN(JJ) 

T=0.0 

DO 104 J=JS,N2 
104 T=C(KK,J1*YCK1,J)+T 

103 YCKl ,JJ)=(Y(K1, JJJ-T)/C(KK,JJ) 

100 Continue 

RETURN 

END 

SUUKUUTINE D3 

C THIS PRUGRAM PERKURMES THE ARTHIMETIC TO DETERMINE WING TIP 

C VELUCITY PUWEK SPECTRUM. SUB... SPEC PERFORMS THE CALCULATION 

C AND FUN...TSPEC EVADURATES ATMOSPHERIC TURBULENCE SPECTRUM. 
C*** N2=NUMbER UF GUST STATIONS 
C*** N=NUMSER FUR DRIVING KkEOUENCIES 
C*** TLsTUkllULENCE LENGTH SCALE 
CUKMUN/TRAN2/RKKC37) 

COMMON /TR AN J/SYC 37, 19,5) 

CuMMUN/EUAT/EEE 
COMPLEX SY 
CALL CUEFl 
CALL COEK2 
HH=19./2. 

DO 101 Jl=l,5 
TL=66.*(24*J1) 

N = 37 

N22=20 

N22=N22«1 
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N2=N22*2 

WK1TEC7 ,202)EEE,TIj,N2 

202 FUKHATUX, ‘Els' ,F10.1, •Tli=* ,F10 .2 , ' N2= ' , 13 ) 

DU 100 1 = 1, 1'J' 
kn=«KK C I ) 

HNU = K^*^1,/BR 

SS=0.0 

U=575. 

Td=T5PEC(SS,RHU,U,TL) 

CaLE SPECCRK,RR»TL,N2,I) 

WRITE (7, 201) RK,RK,TS 
201 FfJRMAT(2X,E14.7,2X,ei4.7,2X,E14.7) 

100 CUiMllNUE 
102 CUivifllJOE 

101 continue 

RETURN 

END 

SUBROUTINE SPEC(RK,RR,TL,N2,1C) 

C THIS SUBROUTINE DETERMINES THE SPECTRUM OF THE WING 
C ’TIP VELOCITY. 

C*** RR=TOTAL AIRPLANE RESPONSE 
C++* RK=REUUCE fREOUENCY riB/U 
C + + * Z=RESPrjNSE TO GUST AT ONE STATION 
C++* N2=NUMBER OF GUST STATIONS 
CuNMohi/TRAn2/RRKC37) 

CUiVMUfJ/TRAN3/SY(37, 19,5) 
dimension YC20.5) ,2140) 

CONILEX Y.CI,2,T,T2,SY 
C++* READ DATA 

C1=CMPLX(0.,1.) 

Bk=19./2. 

022=02/2 

DtL=o6./N2 

U=57b. 

DU 100 J=1,N22 
DO 2 J2=l,5 
Y(J,U2J=SY(1C,J,J2) 

2 continue 

100 CUNTlMiJE 

C + + * UETERMI.jE planes RESPONSE 
DO 101 J=1,N22 
T = 0.0 
T2 = U.O 

Du 102 1=2,5 
T2=Y(J,1)+T2 

IFII .E(J.3.UR.I.EQ.5) T2=T2-2 . *Y ( J , I ) 

102 T=tCJ,l)^T 
ZtJ)=T2+HK+Cl 

101 21U+N22 )=C1+RK*T 

C++* DETERMINE PLNES TOTAL RESONSE 
TT=0.0 

DU 300 1=1, N2 
IS=1-1 
T = 0.0 
Jii=N2-lS 
DO 301 J=1,JN 

301 T=2*REAL(Z(J)*C0NUGC2(J+IS)))+T 
SS=IS*UEL/TL 
RHU=RK + TL/Bi< 
lF(lS.EO.O)T=T/2. 

300 rT=TSPEC (SS , RNU , U , TL) +T+TT 
RR=TT 
RETURN 
END 
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FUlMCTlfUJ TSHEC(5S, RNU,U,Tii) 

c iHiii FurJctiuN uetf;rhines tuubuijEnce cross and power 
c spectrum from the von karmah spectrum function, 

C**» SSsSEP'EKATlOiJ DIVIDED UK Tl. C TURBtiUNCE LENGTH SCALE) 

C*** RfUJ= W*TL/U THE REDUCED FRQUENCT OF TURBULENCE 
C*** U=FL1GHT SPEED UK MEAN WIND SPEED 
C*** TL=TURHULENCE LENGTH SCALE 
IF(SS.EQ.O.O) GU TO SOO 
C CRlJSS SPECTRUM 

Z=SS*SOKT(l .+(1 .339*RNU)**2)/1.339 

TSPEC=TL*.10853/U»(4.78112*SS**(b./3.)/Z4‘*(5,/6. )*BSLUZ) 
+-SS**(ll./3.)/Z**(ll./6.)*BSL2(Z)) 

RETURN 

500 CONTINUE 
C PUrttR SPECTRUM 

Z=(1.339*RNU)**2 

TSPEC=TL*a + C8./3.)*Z)/U+Z)**(ll./6.)/3.14159/U 
RETURN 

Subroutine coefi 

THIS SUliKuUTINE SETS UP THE CUEFFIENETS FOR THE POLYNOMIAL 
APPKGXIHATIOIJ FUR THE HODIFRIED BESSEL FUNCTION OF THE 
SECUfJl' KIND 5/ fa ORDER. 

CuMHui)/Kl 3/A( 10) , B( 10) , A2( 10) 

F=b./fa. 

A( 1 )=1 .0/. 9405612296 
DU 100 1=1,9 

100 A(1+1)=A(I)/1/(F+1) 

F=1 .0-F 

B(l)=l. 0/5. 66756615 
DO 101 1=1,9 

101 B(1 + D=H(I)/I/(F+I-1.0) 

S=4.*(5./6. J**2 

A2(l J=l. 

DO 200 1=1,9 

200 A2(H-n=A2Cl)*CS-(2*l-l)*»2)/8./I 

RETURN 
EN i.) 

FUNCTION bSLKZ) 

THIS FUNCTION EVALURATES THE MODIFRIED BESSEL FUNCTION 
OF THE SECOND KIND 5/6 ORDER 
C0MMUN/K13/A(10) ,B(10) ,A2(10) 

IF(Z.LE.2) GO TO 100 
Y=1 ./Z 

BSLl=SuRT(1.5707*Y)*EXP(-Z)tPOLY(A2,10,Y) 

RETURN 

100 Y=(Z/2.0)<‘*2,0 

R1P=(Z/2.0)**(5,/6.)*POLY(A,10,Y) 

R1N=P0LY(B, 10,Y)/( (Z/2.0)**C5./6.)) 

asLl = ( 3, 14 1/2/SIN (5. 0*3. 141/6,0) )*(RIN-RIP) 

RETURN 

END 

FUNCTION POLYCA.N.Z) 

C THIS FUNCTION DOES THE POLYNOMIAL EVALURATIONS 
DIMENSluN A(N) 

T=A(N)*Z 

NN=N-2 

Du 100 1=1, NH 
100 T=( T-fAlN-I) )*Z 

PULY=T+A(1) 

RETURN 

END 
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FUNCl’lU.M BSL.2CZ3 

C THla FUrJCTIUN tVALURATES THE MODIFRIED BESSEL FUNCTION OF THE 
C SECUUD KIND 11/6 ORDER. 

COMMON /'K2 3/E ( 10),GC10),E2(10) 

lF(Z.Lh.2) GO TO 100 

Y=l./Z 

BSL2=SQRTC1 .5707*Y1*EXP(-Z)*POLY(E2,10,Y) 

RETURN 

100 Y=(Z/2.0)*+2.0 

RIP=CZ/2.0)**Cll./6. )*POLY(E, 10, Y) 

KIK = PUDY(G, 10,Y)/U2./2.0)**(11./6. )) 

BSD2= C 3. 141 /2/SIN (11.0*3, 141/6, 0) )*CRIN-RIP) 

RETURN 
END 

SUBROUTINE COEF2 

THIS SIJHROUTINE SETS UP THE COEFFICIENTS FUR THE POLYNOMIAL 
APPROXIMATIONS OF THE MODIFRIED BESSEL FUNCTION OF THE 
SECOND KIND 11/6 ORDER, 

COKMON/K23/EC10) ,G(10),E2(10) 

F=ll ./tj. 

ONE OVER T[)E GAMMA VALUE OF ItORDER ♦♦*>*i*** 

E( U=1 .0/1 .724362254 
DO 100 1=1,9 

100 E(i+l)=E(l)/I/(F+I) 

(jNE minus the ORDER OF THE MQDFRIED BESSEL ** 

F=1 .0-F 

G(l)=l .O/C-6. 68107938) 

DO 101 1=1,9 

101 G(I + n=G(I)/I/(F+I-1.0) 

S=4*(ll./6.)**2 

E2U) = 1. 

DU 200 1=1,9 . - 

200 E2(I+l)=E2Cl)*(S-C2*I-l)**2)/8,/I 

RETURN 
END 
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APPENDIX B 


ESTIMATION OF WING STIFFNESS 

A major parameter of this study is the beam stiffness. Because 
structural information about the B-57 could not be located, this 
parameter had to be estimated. A static analysis is used in this appen- 
dix to estimate the beam stiffness. We assume that the deflection at 
the wing tip is one foot for a 10 g loading of the aircraft in wartime 
operation (meaning fuel tanks completely full and aircraft loaded with 
bombs both in the fuselage and on the wings). The wing is modeled as 
a cantilever beam. The loading on the wing, including structural weight, 
fuel and bombs, is depicted: 



LOADING OF B-57 WING 
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The values for the loading diagram were determined from the 
Air Force Basic Flight Manual for the B-57. The loading can be written 
in functional form using singularity functions: 

£(y) = -150<y>° - 500<y - 11>° - 550<y - 22>° + 900<y - 25>° 

- 3325<y - 29>° + Q<y>'^ 

where Q represents a hypothetical force at the wing tip and is used for 
Castigl iano' s Theorem to calculate the deflection at the wing tip. The 
loading equation is integrated twice to give the moment equation: 

2M(y) = -150<y>^ - 500<y - 11>^ - 550<y - 22>^ + 900<y - 25>^ 

- 3325<y - 29>^ + Q/2<y>° 

The form of the beam stiffness must also be specified: 


_l , 

0 22 

Distance from wing tip (feet) 


--lOe 
-- e El 

33 


FORM OF El FUNCTION 
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e is the characteristic beam stiffness and must be determined from this 
analysis. The internal strain energy may be written: 


U 

n 




22 

• 0 


id.. 


33 
■ 22 



dy 


where n represents the load factor, assumed to be 10 g. According to 
Castigliano's Theorem, the deflection is: 


A = iU 
h 3Q 


33 


mM 

90 

El 


dy 


1 

e 


22 


33 


dy 


lOe 


'22 


dy 


Substituting the equation for the moment and the derivative of the 
moment with respect to the hypothetical force, the equation becomes: 

r 22 


Ae _ 
h 


•'0 



500 

2 


2 '] 

<y - ll> 

J 


dx 


10 


33 


22 


<y>' 


^<y- 11>2 


^ <y - 22>2 


900 


<y - 25>^ - <y - 29>^J dy 


Performing the integration, the equation is approximately: 


^ 9 X 10® 


Assuming a one foot deflection at the wing tip due to a 10 g loading 
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we get e ~ 9 x 10^. Therefore the beam stiffness distribution due to 
design characteristics is estimated as: 


El = 9 X 

10® 

for 

1 A 

111 

El = 9 X 

10^ 

for 

A 1 

111 
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